rm(list = ls())

## 

library(tidyverse)

##

bt <- read_rds("data/data_federal.rds") %>% 
  mutate(east = ifelse(as.numeric(state_id > 11), 1, 0)) %>% 
  filter(year == 2013) %>% 
  dplyr::select(ags, matches('pop'), treated, east) %>%
  mutate(pop_change_rate = (1 - (pop_dec_09/pop_post_census_2011)) * 100)%>% 
  filter(!pop_dec_09 > 20000)

## Bin population 

bt$pop_09_bin <- cut(bt$pop_dec_09, breaks = seq(0, 20000, 2000),
                     dig.lab = 10)

## Population change by bin

bins <- bt %>% group_by(pop_09_bin) %>% 
  summarise(m = mean(pop_change_rate, na.rm = T)) %>% 
  ungroup()

# Figure 2: Effect of census on population change ----

p1 <- ggplot(bins, aes(pop_09_bin, m)) +
  geom_bar(stat = 'identity', fill = 'grey93', color = 'black',
           width = 0.6) +
  theme_bw() +
  geom_hline(yintercept = 0, linetype = 'solid') +
  geom_vline(xintercept = 5.5, linetype = 'dotted') +
  ylab('Change in population after census\n(percentage points)') +
  xlab('Binned pre-census population (2009)') +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
p1

## East / West comparison

bins <- bt %>% group_by(pop_09_bin, east) %>% 
  summarise(m = mean(pop_change_rate, na.rm = T)) %>% 
  ungroup() %>% 
  mutate(east = ifelse(east == 1, "East Germany", "West Germany"))

## Figure A.3: Census cutoff and changes in population, by East/West 

p1 <- ggplot(bins, aes(pop_09_bin, m)) +
  geom_bar(stat = 'identity', fill = 'grey93', color = 'black',
           width = 0.6) +
  theme_bw() +
  geom_hline(yintercept = 0, linetype = 'solid') +
  geom_vline(xintercept = 5.5, linetype = 'dotted') +
  ylab('Change in population\nafter the census (p.p.)') +
  xlab('Binned pre-census population (2009)') +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  facet_wrap(~east)
p1

